par(mfrow=c(2, 2))
data4_6 <- read.csv("data/习题数据（基于R，EXCEL格式）/csv/E4_6.csv")
ts_x <- ts(data4_6$x)

# 时序图
plot(ts_x, type="o", pch=5, col='#39CBB4')


# ADF检验，平稳性检验
# install.packages("aTSA")
library(aTSA)
adf.test(ts_x, nlag=3)

# 白噪声，纯随机检验（<0.05为非白噪声）
for( k in 1:3) print(Box.test(ts_x, lag=6*k, type="Ljung-Box"))

# 自相关图
acf(ts_x, lag.max=30)
pacf(ts_x, lag.max=30)

# 拟合 AR(1) 模型
fit1 <- arima(ts_x, order=c(1, 0, 0), method="ML")
# 拟合 MA(2) 模型
fit2 <- arima(ts_x, order=c(0, 0, 2), method="ML")

# 模型显著性检验（P大于0.05，模型显著成立）
ts.diag(fit1)
ts.diag(fit2)

# 参数显著性检验
# 粗略的
fit1
fit2
# 精确的 构造t统计量，求P值
t = abs(fit1$coef)/sqrt(diag(fit1$var.coef))
pt(t, length(ts_x)-length(fit1$coef), lower.tail = F)

t = abs(fit2$coef)/sqrt(diag(fit2$var.coef))
pt(t, length(ts_x)-length(fit2$coef), lower.tail = F)



#############
data.frame(AIC(fit1), AIC(fit2), BIC(fit1), BIC(fit2))

install.packages("forecast")
library(forecast)
auto.arima(ts_x)
auto.arima(ts_x, ic="bic")




